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ABSTRACT 

We present the results of 2D simulations of the magnetorotational model of a supernova 
explosion. After the core collapse the core consists of rapidly a rotating proto-neutron 
star and a differentially rotating envelope. The toroidal part of the magnetic energy 
generated by the differential rotation grows linearly with time at the initial stage 
of the evolution of the magnetic field. The linear growth of the toroidal magnetic 
field is terminated by the development of magnetohydrodynamic instability, leading 
to drastic acceleration in the growth of magnetic energy. At the moment when the 
magnetic pressure becomes comparable with the gas pressure at the periphery of the 
proto-neutron star ~ 10 — 15km from the star centre the MHD compression wave 
appears and goes through the envelope of the collapsed iron core. It transforms soon 
to the fast MHD shock and produces a supernova explosion. Our simulations give the 
energy of the explosion 0.6 TO 51 ergs. The amount of the mass ejected by the explosion 
is ~ O.14M0. The implicit numerical method, based on the Lagrangian triangular grid 
of variable structure, was used for the simulations. 
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1 INTRODUCTION 

The problem of explaining core collapse supernova explo- 
sions is one of the long standing and fascinating problems 
, in astrophysics. The most popular core collapse explosion 
mechanisms up to now were the prompt explosion caused 
by the bounce shock and the neutrino driven mechanism. 
The ID spherically symmetrical simulations of the core col- 
lapse supernova do not lead to the explosio n (see for ex- 
ample iBurrows et aljll995l iBuras et al.ll2003l and references 
therein). The 2D and 3D simulations of the neutrino driven 
supernova mechanism do not give supernova explosions ei- 
ther with a sufficient level of confidence. The recently im- 
proved models of the core collapse where the neutrino trans- 
port was simulated by solving the B oltzmann equation do 
not explode either feuras et alJl2003ft . 

The important role of the rotation and the magnetic 
fields for the core coll apse supernova was suggested by 
iBisnovatvi-Koganl jl970T) . The idea of the magnetorotational 
mechanism consists of angular momentum transfer outwards 
using a magnetic field, which twists due to the differential 
rotation. The toroidal component of the magnetic field am- 
plifies with time, it leads an increase in magnetic pressure 
and generation of a supernova shock. 
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The first 2D simulations of the collapse of the ro- 
tating magnetized star w ere presented in the paper of 
iLe Blanck fc Wilson! dl9TCt) . with unrealistically large val- 
ues of the magnetic field. The differential rotation and am- 
plification of the magnetic field resulted in the formation of 
the axial jet. A semi-qualitative pictur e of the magnetoro - 
tational explosions was considered by iMeier et al.l il97Et) . 
The result s of 2D simulat ions of the similar pro blem were 
given in bv lOhnishil <1983h and lSvmbalistvl (Il984l) . Recently, 
simulations of the magnetorotational processes in the col- 
laps ed core have been repre sented by iKotake et alj (l2004tl 
and lYamada fc Sawai (l2004h . where the authors performed 
2D simulations of the magnetorotational processes for a 
strongly magnetized core collapse. They assumed that the 
initial poloidal mag netic field is uniform and parallel to 
the rotation axis. In lKotake et alJ (120041) the initial strong 
toroidal magnetic field was also used as the initial condi- 
tion. In their case the toroidal magnetic field was not equal 
to zero at the rotational axis z, which means a formally 
nonphysical situation with infinite current density at the 
rotational axis. The values of the magnetic fields used by 
them are usual for "magnetars" but too large for the ordi- 
nary core collapse supernova. ID simulations of the success- 
ful magnetorotational supernova explosion for a wide range 
of the initial values o f the magnetic field are represented by 
lArdelian et alj lll979l) . 2D simulations of the magnetorota- 
tional explos ion in a rotating mag netized cloud have been 
described bv lArdelian et all |2000). 
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One of the main problems for the numerical simulation 
of the magnetorotational supernova is the smallness of the 
initial magnetic fields. The ratio of the initial magnetic and 
gravitational energies is in the range of 1(T 6 - IO -10 . This 
value characterizes the stiffness of the set of MHD equa- 
tions. We have a problem with strongly varying character- 
istic timescales. The small one is defined by a very high 
sound speed in the center of the core and the large one is a 
characteristic time of the evolution of the magnetic field. In 
such a situation the explicit numerical methods which are 
widely used in astrophysical hydro-simulations require an 
enormously large number of timesteps and a possible loss 
of accuracy due to the large numerical errors. The implicit 
approach should be used in this case. It is well known that 
implicit schemes are free from the Courant timestep limita- 
tion, which is too strong for such problems. 

At the initial stages of the process the toroidal magnetic 
field is produced by a twisting of the radial component due to 
differential rotation and is growing linearly, B v ~ B r uit. The 
linear growth is terminated when the toroidal component 
becomes so strong that magnetohydrodynamic instability 
starts to develop. This instability leads to poloidal motion 
of the matter, which increases the radial component, which 
in turn strongly amplifies the growth of the toroidal field. 
As a result, both components, toroidal and poloidal, start 
to grow almost exponentially, shortening greatly the time 
between the collapse and magnetorotational explosion. 

The differential rotation and amplification of the 
toroidal component of the magnetic field leads to the forma- 
tion of the MHD compression wave moving from the center 
of the star. This compression wave goes along a steeply de- 
creasing density profile and soon transforms into the MHD 
shock wave. This MHD shock is supported by the flux of the 
matter from the central parts of the star. This matter flux 
works like a piston increasing the velocity and the strength 
of the shock. When the MHD shock reaches the outer bound- 
ary of the collapsed iron core, the kinetic energy of the radial 
motion (supernova explosion energy) is ~ 0.6- 10 51 ergs. The 
amount of mass ejected by the shock is about O.14M0 of the 
total mass of the star, 1.2M0. 



2 EQUATIONS OF STATE AND NEUTRINO 
LOSSES 

For the simulations w e use the equation of state used in 
lArdelian et"afl Jl987al) : 
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where 5R is the gas constant equal to 0.83 ■ 10 8 ^-K, a is the 
constant of the radiation density, P is pressure, p is density, 
and T is temperature. In the expression Po(p) the value p 
was identified with the total mass-energy density. For the 
cold degenerated matter e xpression for Pp(p) i s the approxi- 
matio n of the tables from lBavm et all lll97lf) : lMalone et alJ 
Jl975h . 

In the neighborhood of points p = pt in formula Q 
the function P p (p) wa s smoothed in the same way as in 
lArdelian et"afl Jl987ah . to make continuous the derivative 
dPo/dp : 



Po{p) = 



where 



P 6 [pk-i + £k-i,Pk — £fc] 
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p £ [pk - £,k,Pk + k = 1,5, 



Oh = 0(p) = \ - \ sin \^{P - PkYj , 6 = 0.01p fc . 

The specific energy (per mass unit) was defined ther- 
modynamically as: 

e = £o(p) + §3?T + — + e Fe ( P ,T). 

A p 
The value eo(p) is defined by the relation 
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The term Ep e in equation © is responsible for iron 
dissociation. It is used in the following form: 

Eb,Fe ( T — ToF e 
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It was supposed that in region of the iron dissociation the 
iron amount was about 50% of the mass, Eb,Fe = 8- 10 erg 
is the iron binding energy, A = 56 is the iron atomic weight, 
andrrip = 1.67-10" 24 g is the proton mass, T 0Fe = 0.9-10 10 K, 
TiFe = 1-1 ■ 10 10 K. For the numerical calculations formula 
JSJ has been slightly modified (smoothed): 
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The neutrino losses for Urea processe s are u sed in the 
form, taken from iBisnovatvi-Kogan et alJ (|l976), approxi- 



(7) 



mating the table of Iva nova et alJ ( 196! 

/(P,T)= U1 ^, erg.g-.s-, 
1 + (7.1 • 10- 5 pT)s 



1, T< 7, 

ae(T) = { 664.31 + 51.024(T- 20), 7£T<20, (8) 

664.31, T > 20, 

T = T ■ 10~ 9 . 

The neutrino losses from pair annihilation, photo pro- 
duction, and plasma were also taken into account. These 
types of the neutrino losses have been approximat ed by the 
interpolation formulae from ISchindler et alJ dl987f) : 
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Qtot — Qpair H~ Q photo ~\~ Qpla 



(9) 



The three terms in © can be written in the following general 
form: 



Qd = K(p,a)e -—— r- -7 — -— - — r-. 
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For d = photo, K{p, a) — (p/ pz)oT h ; 
For d = plasm, K(p, a) — (p/pz) 3 ; 
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Here, pz = 2 is the number of nucleons per electron. 
Coefficients c, a u and bj for the differe nt losses are given in 
the Table^from Schindle r et al.l (I1987rl . The general formula 
for the neutrino losses in a nontransparent star has been 
written in the form, used bv lArdelian et aP {2fj 



F(p,T) = (f( P ,T) + Q tot )e- 
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The multiplier e _T ° in formula l|llfl . where r v = S v nl v , 
restricts the neutrino flux for non zero depth to neutrino 
interaction with matter t„. The cross-section for this inter- 
action S v was presented in the form: 



Su = 



1Q -44 T 2 



(0.5965- 10 10 ) 2 ' 
the concentration of nucleons is 

m p 

The characteristic length scale /„, which defines the depth 
for the neutrino absorbtion, was taken to be equal to the 
characteristic length of the density variation as: 

P P 



lu 
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The value l v monotonically decreases when moving to the 
outward boundary; its maximum is in the center, ft approxi- 
mately determines the depth of the neutrinoabsorbing mat- 
ter. The multiplier 1/10 in the expression e _r "To was ap- 
plied because in the degenerate matter of the hot neutron 
star only some of the nucleons with the energy near Fermi 
boundary, approximately 1/10, take part in the neutrino 
processes. 



3 BASIC EQUATIONS 

Consider a set of magnetohydrodynamical equations with 
selfgravitation and infinite conductivity: 
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P— I - ] = H • Vv, A$ = 4ttGp, (13) 



p^ + PV ■ v + P F(p,T) =0, 

P = P{p,T),e = e{p,T), 

where i = + v ■ V is the total time derivative, x = 
(r, tp, z), v = (v r , v v ,v z ) is the velocity vector, p is the den- 
sity, P is the pressure, H = (H r , H v , H z ) is the magnetic 
field vector, <E> is the gravitational potential, e is the internal 
energy, G is gravitational constant, H Cg> H is the tensor of 
rank 2, and F(p, T) is the rate of neutrino losses. 

r, ip, and z are spatial Lagrangian coordinates, i.e. 
r = r(r ,(p , and z ,t), <p = (p(r , <fio, Zo, t), and z = 
z(ro, <po, zo, t), where ro, (fo, zo are the initial coordinates of 
material points of the matter. 

Taking into account symmetry assumptions (^- = 0), 
the divergency of the tensor H ® H can be presented in the 
following form: 
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Axial symmetry (^- = 0, r ^ 0) and symmetry to the 
equatorial plane (z ^ 0) are assumed. The problem is solved 
in the restricted domain. At t = the domain is restricted 
by the rotational axis r ^ 0, equatorial plane z ^ 0, and 
outer boundary of the star where the density of the matter 
is zero, while poloidal components of the magnetic field H r , 
and Hz can be non-zero. 

We assume axial and equatorial symmetry (r ^ 0, z ^ 
0). At the rotational axis (r = 0) the following boundary 
conditions are defined: (V$) r = 0, v r = 0. At the equato- 
rial plane (z = 0) the boundary conditions are: (V<J>) 2 = 
0, Vz — 0. At the outer boundary (boundary with vacuum) 
the following condition is defined: P ou t er boundary ~ 

We avoid explicit calculations of the function £o(p) in 
(0J|, because this term is eliminated from 113H due to adia- 
batic equality: 
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determining the fully degenerate part of EOS. Therefore, 
defining 
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The equation for the internal energy in l|130 can be written 
in the following form: 



P^+P*V-v + pF(p,T) =0. 



(15) 



4 DIMENSIONLESS FORM OF EQUATIONS 

For the simulations we rewrite the set of equations 1131 in 
the dimensionless form. The basic scale values are: 



r = 1.35 ■ 10 s cm, p = 10 9 g/cm 3 , 
G = 6.67- 10~ 8 cm 3 /(g-s 2 ). 



(16) 



4 N.V.Ardeljan, G.S.Bisnovatyi-Kogan, and S.G.Moiseenko 



Table 1. 


CoefBcients for formula 1101 fromlSchindler et alJ 
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2.146(-7) 


1.745(20) 
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The dimensional functions can be represented in the follow- 
ing way (values with " tilde sign are dimensionless values): 

r = fr , z = zr , p = ppo, v = vvo, 
t = tt , v r = v r v , v v = v v vq, v z = v z vo, p = ppo, (17) 
T = TTq , $ = l-$ = $47rGp ro, e = ee , H = HH , (18) 
where 

v = y/4nGp r 2 = 3.908 • 10 9 cm/s, 

, r 2 v 2 

to = — , Po = Po^o, ^0 = 75-, 
Vq K 

$0 = 4nGp ro, e = Vq, H = y^o" = x to 1 pa~ 1 - 

Taking into account 1151 the set of basic equations I13H can 
be written in the following dimensionless form (the tilde sign 
being omitted here): 
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5 INITIAL CONDITIONS 



We have used the model I from lArdelian et alJ (^87^) as ini- 
tial conditions . At first we calculated spherically symmetri- 
cal stationary model with central density p c = 4.5-10 9 g/cm 3 . 
The value of the central density corresponds to the maxi- 
mum in the dependence of the stellar mass on the central 
density Ms(p c ) for T = 0. The mass of such a spherical star 
is M = 1,OO42M . 

To define the initial model the density (and hence the 
mass) of the star was increased by 20% at every point. 
The temperature in the star was defined by the relation: 
T = 8p 2/3 , where 5 = 1(K • cm 2 • g~ 2/3 ). At the initial time 
moment t = 0, the rigid-body rotation with the angular ve- 
locity u) = 2.519 c _1 (rotational period is r = 2.496 c) was 
accepted. At t = we suppose also that v r — v z = 0. The 



initial rotational energy is 0.571% of the gravitational en- 
ergy, and the initial internal energy including the energy of 
degeneracy is 72.7% of the gravitational energy. The defined 
initial model is unstable against collapse, and immediately 
after the beginning of the calculations it starts to contract. 

The initial magnetic field (the initial magnetic field 
which was "turned on" after the collapse stage) was de fined 
exactly in the same way as in lArdelian et all tOOd) . We 
define the toroidal current j v in the following form: 



ft forz ^ 0, r 2 + z 2 < 0.025 2 , 



3<p = < 3% for 2^0, 
I for 



r 2 +z 2 ^ 0.025 2 , 



(20) 



where 



i u — A ■ 
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Aj is a coefficient used for adjusting the values of the ini- 
tial toroidal current and, hence, the magnetic field. After 
the collapse of the core and obtaining the differentially ro- 
tating stationary configuration without a magnetic field, we 
calculated the initial poloidal magnetic field H r o,HzO us- 
ing the Bio-Savara law. The calculated magnetic field is 
divergence-free but not force-free and not balanced with 
other fo rces. To make it ba lanced we use the following 
method jArdelian et alj|200(j) : we "turn on" the poloidal 
magnetic field H t q, H z q, but "switch off" the equation for 
the evolution of the toroidal component H v in l|13|l . In fact 
this means that we define H v = 0, —^f- = 0. From the 
physical point of view it means that we allow magnetic field 
lines to slip through the matter of the star in the toroidal 
direction <p. After "turning on" such a field, we let the cloud 
come to a steady state, where magnetic forces connected 
with the purely poloidal field are balanced by other forces. 
The calculated balanced configuration has a magnetic field 
of quadrupole-like symmetry. The relation of the magnetic 
energy of the star to its gravitational energy at the moment 
of the formation of the balanced poloidal magnetic field is 
10" 6 . 
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6 NUMERICAL METHOD 

The implicit operator-difference completely conservative 
scheme on the triangular grid of variable structure was used 
for the numerical modeling of the problem of the magnetoro- 
tational supern ova explosion. The sche me was suggested and 
investigated bv lArdelian et all il987bf) and in earlier papers 
by these authors. 

The solution of the problem is a sequence of time steps. 
Calculation of every time step can be divided into two parts. 

The first part is the calculation of the values of the func- 
tions on the next time level using the implicit completely 
conservative operator-differenc e scheme on the triangu- 
lar gri d in Lagrangian v ariables lArdelian fc Kosmachevskiil 
(Il995l) . lArdelian et al ] il987tJ> . The coordinates of grid 
knots are changing at this stage. 

The second part is an analysis of the quality of the 
grid, its improvement and adaptation (grid reconstruction). 
The improvement of the quality of the grid is necessary be- 
cause of the appearance of "poor" cells, i.e. triangles which 
strongly deviate from equilateral triangles. The dynamical 
adaptation of the grid allows us to concentrate the grid in 
the regions of the computational domain where spatial reso- 
lution needs to be increased and rarefy the grid in the regions 
where the flow is smooth. It allows us to reduce significantly 
the dimensionality of the grid and hence strongly reduce the 
computation time. 

The grid reconstruction procedure itself consists of the 
following stages. The first stage is a local correction of the 
structure of the grid. The second stage is a calculation of 
the values of the functions defined in cells in the regions of 
the corrected structure. 

The local correction of the structure of the grid can be 
d one using the fo l lowing three local operations (see details 
in lArdelian et alJ Jl996f) 1: 

(i) replacement of the diagonal of the quadrangle formed 
by two triangles by another diagonal; 

(ii) joining up of two neighboring grid knots; 

(iii) adding of the knot at the middle of the cell side which 
connects two knots. 

The improvement of the grid structure is made using the 
first two operations. The grid adaptation is made by appli- 
cation of the local operations (II) and (III). The plots which 
clearly explain these operations are given in lArdelian et ail 
(1996). 

The values of functions defined in cells and knots in- 
volved in grid structure modification are calculated at the 
second stage of the grid reconstruction. 

The application of simple interpolation for the calcula- 
tions of new values of the grid function leads to the violation 
of the conservation laws and adds significant errors in the 
regions with high gradients (for example at shock waves). 
The goal of the calculation of the new values of functions 
is to minimize numerical errors introduced by this proce- 
dure. To achieve this goal not only the error in the values 
of the functions but also on their gradients have to be min- 
imized. It is important to fulfill conservation laws (mass, 
momentum, energy, magnetic flux) in the vicinity of the lo- 
cal grid reconstruction. The method for the calculation of 
these new values of grid functions is based on a minimiza- 
tion of the functionals containing the values of the func- 



tions, its gradients, and grid analo gs of the conservation laws 
(lArdelian fe Kosmachevskiilll995h . 

For the solution we have used the method of the con- 
ditional minimization of the functionals guaranteeing exact 
fulfilment of the conservation laws. 

It is important to fulfill conservation laws for the solu- 
tion of the collapse and magnetorotational supernova explo- 
sion problems, because a large number of the time steps need 
to be made. In such a situation even slight violation of the 
conservation laws at a time step could lead to a significant 
growth in the errors and hence to a qualitative distortion of 
the results. 

The method of calculating hydrodynamical values 
p and P during grid recons truction is taken from 
lArdelian fc Kosmachevskiil (^)95) . To obtain a better preci- 
sion some changes have been made in the calculation of the 
magnetic field components in comparison w ith our calcula- 
tions of the collapse of magnetized cloud llArdelian et alJ 
(2000)). Earlier, for the calculation of the new values of 
the magnetic field components we conserved the sum of the 
products of the magnetic field components values and the 
volumes of the cells in the vicinity of the reconstructed part 
of the grid. In this paper we conserve the poloidal magnetic 
energy and the toroidal magnetic flux in the vicinity of the 
reconstructed part of the grid. 

The grid reconstruction procedure allows us not only 
to "correct" the Lagrangian grid, but also to dynamically 
adapt it using different criteria for the grid in different parts 
of the computational domain. The grid can be refined in the 
regions where it is necessary and therefore we can increase 
the accuracy of the calculations. It is possible to rarefy the 
grid in those parts of the computational domain where the 
flow is "smooth". The procedure of the rarefying of the grid 
allows us to significantly reduce the dimension of the grid 
while preserving the same accuracy for the numerical solu- 
tion. 

Geometrical criteria can be used for the grid adapta- 
tion (i.e. restrictions on the length of the cell size, which 
are defined by the cell coordinates only), but such adap- 
tation criteria are suitable for the flows of simple or easily 
predictable structures only. 

For grid adaptation in the case of complicated flows or 
flows with unknown structures it is better to apply dynami- 
cal criteria which are defined by the solution behavior. The 
dynamical cr iteria for the gr i d ada ptation applied here was 
suggested bv lArdelian et al] ]l996h . 

The characteristic length of the cell side It was chosen as 
a local criterion for the grid reconstruction. As an example, 
consider the criterion where 1^ is a function of p and gradp. 
Let's introduce the function: 



/(p, gradp) = 



(p + e) 71 (|gradp|+e)72 



(21) 



where < e << 1, a ^ 0, /3 > 0, a + f3 = 1. 71, 72 are 
power indexes, and gradp is the grid analog of the density 
gradient. In limited cases: 

a = 1, P — 0, / depends on the density only; 

a = 0, /3 = 1, / depends on the gradient of the density 

only. 

Let be the total number of the grid cells. The char- 
acteristic length of the side Ik of the cell with the number k 
will be calculated as a function of the density, the gradient 
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of the density, and the coordinates 
following formula 

f(Pk, gradp fe ) 



Ik = 2 



Sk = 



(implicitly) by the 



(22) 



J2 f(Pn, gradp„ 



Here, S is a square of the computational domain, which 
consists of the triangular cells. The summation in the de- 

N 

nominator is made for all grid cells. Note that S = Sh, 

k = l 

where Sk is equal to the square of the equilateral trian- 
gle with the length of the side Ik- In our calculations the 
function / in formula 1211 was used in the following form 
f(pk) = l/( P k+e) - 5 . 

The criterion described was applied in the following 
way. In the case where the length of the side of the cell i 
is larger than 2lk, a new knot is added to the middle of this 
side of the cell. In the case where the length of the side of 
the cell i is less than Q.Hk, the operation of the joining up 
of these knots is applied. The application of the dynamical 
adaptation criterion described allowed us not only to adapt 
the grid to the specialities of the solution but also to pro- 
vide an acceptable accuracy of the calculations with a small 
fluctuation of the total number of grid knots and cells. In 
the calculations described the total number of knots (5000) 
and cells (10000) was violated not more than 5%. 

At the moment of the maximum compression (for the 
collapse problem) the minimal size of the cell side is so small 
that application of the uniform grid with the same spatial 
resolution would require a grid with the dimension ~ 1000 x 
1000 (!) cells. 

The calculation of the grav itational potential in this pa- 
per as in lArdelian et alJ 1120041) in the frames of the applied 
numerical method is ma de on the base of the f inite e lement 
method of higher order JZenkevich fc Morgan! Jl983t) ). This 
procedure allowed us to increase the accuracy of the calcu- 
lations and to eliminate the loss of approximation near the 
z axis. 

T he linear artificial viscosity (e.e.. ISamarskii fc Popovl 
1992) was introduced into the numerical scheme for the 
shock capturing. 



7 RESULTS 

The simulation of the magnetorotational supernova is di- 
vided into two separate parts. The first part is the core col- 
lapse simulation, and the second is the "switching on" of 
the initial poloidal magnetic field following amplification of 
its toroidal component H v , which is finished by the magne- 
torotational explosion. The duration of the second phase is 
defined by the value of the initial magnetic field. The weaker 
the initial poloidal magnetic field, the longer the stage of the 
amplification of the toroidal component until the supernova 
explosion. The results of ID simu lations of the magnetoro- 
tati onal mechanism made bv lBisnovatvi-Kogan et all jl976t) 
and lArdelian et all il979fl show that the time from the be- 
ginning of the magnetic field evolution to the explosion is 

proportional to 1/ \ 




Figure 1. Time evolution of the neutrino luminosity 

Moore 

J f{p> T)dm during collapse. 




7.1 Core collapse simulation 

The core collapse sim ulation was described in detail in 
lArdelian et alJ i2004f) . In the present paper we describe 
briefly the obtained results. The ratios between the initial 
rotational and gravitational energies and between the inter- 
nal and gravitational energies of the star are the following: 

' ' ' = 0.0057, = 0.727. 



E 



E a 



where E mag0 and E grav0 are 
initial magnetic and gravitational energies of the star, re- 
spectively. 



Soon after the beginning of the contraction at t — 0.1377s 
at the distance 6 • 10 5 cm the bounce shock wave appears. 
Behind the shock front the temperature of the matter rises 
sharply and neutrino losses are "switching on". In Fig.^the 
time evolution of the neutrino losses is presented. The mat- 
ter of the star from the outside of the shock wave continues 
to collapse towards the center of the star. 

At t — 0.1424s the density in the center of the star 
reaches its maximum value p c ,max = 5.655 • 10 14 g/cm 3 . The 
matter of the envelope which has passed through the bounce 
shock wave forms the core (proto-neutron star). Behind the 
shock front the intensive mixing of the matter takes place. 

The shock wave moves through the envelope and at 
t = 0.2565s reaches the outer boundary of our computational 
domain near rotational axis z. The shock leads to the ejec- 
tion of 0.041% of the core mass and 0.0012% (2.960- 10 48 erg) 
of the gravitational energy of the star. The particles of the 
matter are treated as "ejected" when their kinetic energy 
becomes larger than their gravitational energy and the ve- 
locity vector is directed from the star center (r = 0, z = 0). 
The amounts of the ejected mass and energy are too small 
to explain the c ore collapse sup e rnova explosion. It should 
be noted that in lJanka fc Plewal (120021) . where the neutrino 
losses were calculated using the solution of the Boltzmann 
equation, the shock wave does not lead to the ejection of the 
matter. At the final stage of the core collapse (t = 0.261s.) 
we obtain a differentially rotating configuration. The central 
proto-neutron star with a radius ~ 12.8km rotates almost 
rigidly with the rotational period 0.00152s. The angular ve- 
locity rapidly decreases with the increase in the distance 
from the star center. The particles of the matter situated at 
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R [10*cni] 

Figure 2. Dependence of the rotational period (in s.) from r 
coordinate in the equatorial plane (z = 0) after the collapse for 
t = 0.261s. 

the outer boundary in the equatorial plain rotate with the 
period ~ 35s (Fig.EJ. 

7.2 Magnetorotational explosion 

The results of the collapse simulations show that the 
amounts of the mass and energy ejected by the bounce shock 
wave are too small to explain the supernova explosion. 

The small ejecting part of the envelope prevents us from 
continuing the calculations because of significant increase in 
the square of the computational domain. To overcome this 
problem the outer part of the envelope was stopped by arti- 
ficially reducing its velocity (kinetic energy). This artificial 
method does not influence the central parts of the collapsed 
core. 

After "stopping" the ejecting part of the envelope we 
run our code without the magnetic field to be sure that the 
configuration formed is not changing in time significantly. 
We reach the stage where the poloidal kinetic energy of the 
star is less than 0.00001% of the gravitational energy of the 
star and remains below that value for 1000 time steps. The 
shape of the envelope of the star does not change signifi- 
cantly during this test. 

At that stage we "switch on" the initial poloidal mag- 
netic field, defined by the current l)20jl . The energy of this 
magnetic field was taken to be equal to ~ 10~ 6 of the grav- 
itational energy of the star at the moment of including the 
magnetic field in our simulations. The poloidal magnetic 
field, defined by the current l)20|l . is divergence- free, but is 
not force- free, and "switching on" that field can lead to the 
artificial violation of the equilibrium of the star. To reach 
a steady state of differentially rotating configuration with 
the balanced p oloidal magneti c field we exactly follow the 
procedure from lArdelian et al.l l)2000l) . 

The balanced configuration calculated has the magnetic 
field of quadrupole-like symmetry. After the formation of the 
balanced poloidal magnetized configuration we "switch on" 
the equation for the evolution of the toroidal component H v . 

We calculate the evolution of the magnetic field af- 
ter the collapse stage, because at the developed stage the 




Figure 3. Time evolution of the toroidal and poloidal parts of 
the magnetic energy during magnetorotational explosion. 



collapse is rather short in time and newly forming proto- 
neutron star rotates not as differentially as at the end of the 
collapse. During the collapse the forming proto-neutron star 
makes only a few revolutions and the magnetic field, which 
is initially weak, does not have a significant influence on the 
flow in the star. 

At the moment of "switching on" the toroidal magnetic 
field we start a counting the time anew. 

At the beginning of the simulations the toroidal com- 
ponent of the magnetic field grows linearly with the time at 
the periphery of the proto-neutron star. The energy of the 
toroidal magnetic field grows as a quadratic function (Fig. 
0. At the developed stage of the H v evolution (t = 0.04s) 
the poloidal magnetic energy begins to grow much faster due 
to de veloping magnetorotational instability dAkivama et alJ 
2003) leading also to a rapid growth of the poloidal com- 
ponents. Due to the quadrupole-like type of the symme- 
try of the initial magnetic field the generated toroidal field 
H v has two extremes. The first is at the equatorial plane 
and the second is out of the equatorial plane at the pe- 
riphery of the proto-neutron star closer to the axis of ro- 
tation z. The two extremes have different signs of H v . In 
Fig-Elthe distribution of the H v is plotted at 0.012s. These 
extremes approximate ly correspond to the extremes of the 
term rH ■ grad{V v /r) iArdelian et al.ll2000l) in the equation 
for the evolution of the H v , because the star is in a steady 
state condition, and only this term in this case determines 
the evolution of H v . 

The maximal value of the H v reached during the ampli- 
fication of the toroidal field phase is ~ 2.5- 10 16 G. This maxi- 
mum is situated at the equatorial plane at a distance 1.7- 10 6 
cm from the center of the star, and it is reached at t — 0.058s. 
The toroidal part of the magnetic energy decreases with time 
after reaching its maximal value of 4.8- 10 50 ergs at t = 0.12s. 
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Figure 4. The toroidal field H v distribution after 0.012s of the 
beginning of the evolution of the H v in the central part of the 
star. 



The poloidal magnetic energy at the developed explosion 
stage is w 2.5 • 10 50 erg and keeps this value until the end of 
our simulations (Fig. . 

The magnetic pressure is highest in the regions where 
the H^p reaches its extremal values. At these regions the \Hip\ 
is approximately 100 times higher than the absolute value 
of the poloidal field {VH? + H%). 

The angular momentum is subtracted by the mag- 
netic torque from the proto-neutron star. The magnetic field 
works as a "transmission belt" for the angular momentum. 
The envelope of the star starts to blow up slowly. The con- 
traction wave appears at the periphery of the proto-neutron 
star near the extremes of the H v . This contraction wave 
moves along a steeply decreasing density profile. The am- 
plitude of this wave rapidly grows, and after a short time 
it transforms into the fast MHD shock wave. The growing 
toroidal magnetic field due to the differential rotation works 
as a piston for the originated MHD shock. Time evolution of 
the velocity field v T ,v z , specific angular momentum v v r, and 
temperature T for the time moments t = 0.07, 0.20, 0.30s is 
given in Fig. |SJ and Fig. HJ 

The flow behind the initial MHD shock is very inhomo- 
geneous. Possible reasons of this inhomogeneity are: 

• the aftershock behavior of the gas flow in the presence 
of the gravitational field, 

• the magnetorotational instability appearing at the pe- 
riphery of the proto-neutron star. 

The oscillating structure of the flow behind the shock 
wave in the gravitational field (without magnetic fields) was 
investigated an alytically in a ID case in the acoustical ap- 
proximation bv iLambl lll909l) and numerically for ID gravi- 
tational gas dynamics bv lKosovichev fc Popov! lll979T) . Sim- 
ilar oscillating structure was found in our 2D simulations of 
the collapse of the ra pidly rotating cold protostellar cloud 
iArdelian et aljll996t) . The main reason for this effect is a 
change of dispersion properties of the matter in the presence 
of a gravitational field. 

During the evolution of the magnetorotational explo- 
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Figure 5. Time evolution of the specific angular momentum v v r 
for the time moments t = 0.07s, 0.20s, 0.30s. The darker parts of 
the plots correspond to the higher specific angular momentum. 
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Figure 7. Time dependence of the neutrino luminosity 

Mcore 

J f(p, T)dm during the magnetorotational explosion. 




sion the proto- neutron star loses angular momentum, con- 
tracts and slightly increases its rotation. This leads to the 
further amplification of the toroidal magnetic field in this 
region. In other words it means that the "piston" continues 
to push the matter from the central object to the envelope 
and to support the MHD shock. The continuous support of 
the supernova shock in the magnetorotational mechanism is 
the main qualitative difference from the prompt shock and 
neutrino driven supernova mechanisms. 

Due to the quadrupole-like type of the symmetry of 
the initial magnetic field the MHD shock is stronger, and it 
moves faster near the equatorial plane z = 0. The matter 
of the envelope of the star is ejected preferably near the 
equatorial plane. 

The formation of the MHD shock and its propagation 
leads to the secondary neutrino losses jump (Fig. |7J, but 
now it is much weaker than it was at the collapse stage. 

The integral neutrino losses as a function of time (dur- 
ing magnetorotational stage) is given in the Fig. |S] 

Time dependence of the ejected mass of the star is given 
in the Fig. [5] In Fig. HOI the time dependence of the ejected 
energy is represented. The particle is considered as ejected if 
its kinetic energy is greater than its potential energy and its 
velocity vector is directed from the center. The results of our 
simulations show that during magnetorotational explosion 
~ 0.14 Mq of the mass and ~ 0.6 • 10 ergs are ejected. 

The MHD sho ck which produces th e supernova is a fast 
MHD shock (as in lArdelian et al J [20001. because its veloc- 
ity is larger than fast magnetic sound speed in the upstream 
flow, while the velocity of the slow MHD shock in the down- 
stream flow is between Alfvenic and slow magnetic sound 
speeds. 

The magnetic field transforms part (~ 10%) of the rota- 
tional energy of the star to the radial kinetic energy (explo- 
sion energy). The time dependence of the rotational energy 
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Figure 8. The integral neutrino losses as a function of time dur- 
ing the magnetorotational explosion. 
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Figure 9. Time dependence of the ejected mass in relation to 
Mq during the magnetorotational explosion. 



and the poloidal part of the kinetic energy are given in Fig. 

EH 

During the magnetorotational explosion the star loses 
a significant part of its rotational energy. The rotational 
energy of the star is transformed not only into the explosion 
(kinetic energy of radial motion) but is lost in the form of 
neutrino emission and partially changes the total energy of 
the star (Fig. 1121 . The core rotates more slowly now, and 
it leads to a deeper contraction and some heating of the 
proto-neutron star. 

We stop calculations at t ~ 9.45s. At this stage the 
proto-neutron star rotates with the period ~ 0.006s. The 
absolute value of the poloidal magnetic field at the periphery 
of the proto-neutron star (at the equatorial plane, ~ 10km 
from the center of the star) is VH'i+ H'i « 2 • 10 14 G. 
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Figure 10. Time dependence of the ejected energy during mag- 
netorotational exDlosion. 
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Figure 11. Time dependence of the rotational energy and 
poloidal kinetic energy of the star during magnetorotational ex- 
plosion. 



8 APPEARANCE OF THE 

MAGNETOROTATIONAL INSTABILITY IN 
2-D PICTURE 

Magnetorotational instability (MRI) in the magneti z ed sta r 
with differential rotation was analyzed by ISpruitJ J2002I) . 
who indicated to dynamo action, accompanies a develop- 
ment of such instability. In our axial ly symmetric 2 D simu- 
lations dynamo action is prohibited Jcowlinelll957r) . never- 
theless the development of the MRI takes place. It was in- 
ve stigated in 2D numerical ly fo r the case of accretion discs 
bv lHawlev fc Balbusl <199lft and lFromang et alj J2004 . The 
possibility of the appearanc e of su ch instability was also 
mentioned by IColeate et al.l l)l99ft) in application to the 
SN 1987a. The qualitative picture of the MRI in 2D is the 
following. At the first stages the differential rotation leads to 
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Figure 12. Time dependence of the gravitational and internal 
energies of the star during magnetorotational explosion. 





Figure 13. Qualitative picture of the development of the mag- 
netorotational instability. 



a linear growth of the toroidal field described qualitatively 

as 



dt 



dr J 



(23) 



The right side is constant at the initial stage of the process. 
When toroidal field reaches its critical value H v *, the MRI 
instability starts to develop. As follows from our calcula- 
tions the critical value corresponds to the relation between 
toroidal magnetic and internal energy densities as 



H* 



-Ep. 



(24) 



8tt 10 

Appearance of MRI is characterized by formation of multi- 
ple poloidal differentially rotating vortexes, which twist the 
initial poloidal field leading to its amplification according to 



dH r 
~dT 



dw v 
~df 



(25) 



where I is the coordinate, directed along the vortex radius, 
u v is the angular velocity of the poloidal vortex. Qualita- 
tively the poloidal field amplification due to the vortexes in- 
duced by MRI is shown in the Fig. 1131 The enhanced poloidal 
field immediately starts to take part in the toroidal field am- 
plification according to 1231 1. With further growing of H v the 
poloidal vortex speed increases. Our calculations give the 
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values of uj v = 0.0132s" 1 at \H V \ = 2.46- 10 15 G correspond- 
ing to t = 0.041 s and lj v = 0.052 s _1 at \H V \ = 4.25 • 10 15 G 
corresponding to t — 0.052s" 1 for the same Lagrangian par- 
ticle. In general we may approximate the value in brackets 
12511 by linear function on the value (H v — _ff*) as 

(^l)=a(H v -H;). (26) 

H r (or H po i oi dai in general case) in 12311 is not constant any- 
more after onset of MRI, but described by 125H and 126H . 
Assuming for the simplicity that (r^r) = A is constant 
during the first stages of MRI, and taking H* as constant 
we come to the following equation: 

^ (H v -H;) = AH p0 a(H v - H* v ), (27) 
giving the exponential growth of magnetic field components 



H V = H;+ H p0 eV AaH p° <*-**>, 
Hr = Hr0 + H ^J^ 2 (e^ (t - n - l) , 

taking H p o as seed field for development of MRI. 

The above toy model shows the way of development of 
MRI instability in 2D case. In this case there is no direct in- 
fluence of toroidal field onto a poloidal one, therefore there is 
no dynamo action. Nevertheless, both components grow ex- 
ponentially, chaotically for the poloidal field and both chaot- 
ically and regularly for the toroidal one. Magnetorotational 
explosion is produced by both, regular and chaotic magnetic 
fields. After rapid neutrino cooling and damping of the hy- 
drodynamic motion, the chaotic component remains frozen 
into the matter due to high electric conductivity, and chaotic 
magnetic field could be the source of magnetar activit y in 
soft gamma repeaters iThompson fc Beloborodovl i2004T) . 



9 DISCUSSION 

Our calculations have shown that the magnetorotational 
mechanism produces enough energy (0.6 • 10 51 ergs) to ex- 
plain the core-collapse supernova of types II and lb. The Ic 
type is probably more energetic and could be connected with 
the collapse of much more massive cores of a few tens of solar 
masses. The value of the energy production in this particu- 
lar variant may be considered as a basic one, which may be 
several times larger or smaller, depending on the mass of the 
collapsing core, magnetic field intensity and configuration, 
and rotational energy of the presupernova. 

The magnetorotational explosion of a supernova star is 
divided into three stages: 

(i) linear growth of the toroidal magnetic field due to the 
twisting of magnetic field lines, 

(ii) exponential growth of the toroidal and poloidal mag- 
netic fields due t o the development of magnetohydrody - 
namic instabilities l lDunge\ll95klTavleJl97alSpruiil200al . 

(iii) formation of the MHD shock wave and magnetoro- 
tational explosion. 

The resulting newborn neutron star is characterized by 
relatively slow rotation: 6 ms in comparison to sub-ms values 



of the critical rotation. The resulting period could grow with 
increasing efficiency of the explosion. 

The toroidal magnetic energy of the young neutron star 
decreases more rapidly than the poloidal component which 
tends to the constant value (see Fig.[3J. Possibility of the dy- 
namo action in a similar situation, leading to an increase in 
the large scal e poloidal as w ell as toroidal components, was 
considered bv lSpruitl i2002f) . In our calcul ations the dyna mo 
action is prohibited due to 2D geometry jCowlinalll957l) . 

Different variants of magnetorotational instability in as- 
troph ysics were investigated in connection w ith jet forma- 
tion jLovelacelll976l fMeier fc Nakamurall2003j) . and genera- 
tion of turbulence in the accretion disks fealbus fc Hawlevl 
1998). In the last case the instability of the uniform field, 
threading the differentially rotating disk is co nsidered, which 
was investigated ear lier in plasma physics <IVelikhovlll959] . 
Chandrasckhar 1981). 

The large values of the remaining magnetic field in 
our case are related to the chaotic, nonregular components, 
which have zeroaverage magnetic flux and could disappear 
by field annihilation even at large conductivity. The exis- 
tence of a large chaotic component of magnetic field on a 
neutron star may last long due to very high conductivity. 
It may be related the magnetars model or soft gamma re- 
peaters, in which the radiated energy comes from the mag- 
netic field annihilation. Inside the neutron star the regular 
toroidal, dipole, or quadrupole poloidal components of the 
magnetic field could remain for a long time, much longer 
than their chaotic values. 
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